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ABSTRACT 

We describe an accurate, one-dimensional, spherically symmetric, Lagrangian 
hydrodynamics/gravity code, designed to study the effects of radiative cooling and 
photo-ionization on the formation of protogalaxies. The code can treat an arbitrary 
number of fluid shells (representing baryons) and collisionless shells (representing cold 
dark matter). As a test of the code, we reproduce analytic solutions for the pulsation 
behavior of a polytrope and for the self-similar collapse of a spherically symmetric, 
cosmological perturbation. In this paper, we concentrate on the effects of radiative 
cooling, examining the ability of collapsing perturbations to cool within the age of the 
universe. In contrast to some studies based on order-of-magnitude estimates, we find 
that cooling arguments alone cannot explain the sharp upper cutoff observed in the 
galaxy luminosity function. 

Subject headings: Galaxies:formation, Hydrodynamics, Methods: numerical 

1. Introduction 

The leading cosmological theories imply that galaxies form by the collapse of primordial 
density fluctuations. The gravitational evolution of collisionless matter can be followed by various 
dynamical approximations, or, in the strongly non-linear regime, by N-body simulations. However, 
gas dynamical effects such as shocks and radiative cooling must play an essential role in the 
formation of galaxies, since gas must cool and condense inside dark matter halos before it can 
form stars. 

There have been two quite different approaches to this theoretical problem. One, going back 
to the pioneering work of Binney (1977), Silk (1977), and Rees & Ostriker (1977) (hereafter 
RO), uses simple analytic estimates: typically, one computes the characteristic density and virial 
temperature of a dark halo assuming a spherical collapse model, then asks whether gas at this 
density and temperature can cool within a dynamical time, or within a Hubble time. Combined 
with extended versions of the Press-Schechter (1974) formalism, these methods can yield detailed 
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predictions for properties and evolution of the galaxy population (White & Rees 1978; White 
& Frenk 1991; Kauffman, White & Guiderdoni 1993; Cole et al. 1994). The second approach, 
which has become computationally practical only within the last few years, is to incorporate gas 
dynamics directly into three-dimensional numerical simulations (e.g. Katz & Gunn 1991; Cen & 
Ostriker 1992, 1993; Katz, Hernquist & Weinberg 1992; Evrard, Summers & Davis 1994; Steinmetz 
& Muller 1994). 

In this paper and those that follow it, we will take an intermediate path, modeling the 
collapse of individual perturbations with a one-dimensional, Lagrangian, gravity/hydro code. The 
code evolves a mixture of gas and collisionless dark matter, elements of which are represented by 
concentric, spherical shells. The gas responds to gravity and pressure forces; it can be heated by 
adiabatic compression, by shocks, and by energy input from a photo-ionizing background, and it 
can cool by a variety of atomic radiative processes. The collisionless dark matter responds only 
to gravitational forces. While we focus in this paper on galaxy-scale collapses assuming spherical 
symmetry, the code is also well suited to studies of Lyman-alpha clouds, and it can easily be 
adapted to planar or cylindrical symmetry. 

Larson (1969,1974) studied galaxy formation with spherically symmetric simulations more 
than 20 years ago. However, the intervening years have seen great changes in the theoretical 
underpinnings of galaxy formation - especially the introduction of dark matter and the 
development of physically motivated initial conditions - and they have seen great improvements 
in computational algorithms and hardware, so there is plenty of reason to revisit this approach. 
Our calculations include radiative cooling in the gas component and gravitational interactions 
with a collisionless component, and we adopt initial conditions appropriate to Gaussian random 
fluctuations, as might be produced by inflation in the early universe. Instead of Larson's 
Eulerian-grid approach, we adopt a Lagrangian representation of the gas and dark matter, 
which provides much higher spatial resolution in the central, high-density regions of a collapsing 
protogalaxy. The manyfold increase in computer power allows us to use large numbers of fluid 
elements and to perform faster searches in parameter space. 

Our numerical approach — one-dimensional, Lagrangian hydrodynamics with radiative 
cooling — is similar to that used by Thomas (1988) in his models of cooling flow galaxies and 
by Shapiro & Struck-Marcell (1985) in their studies of "pancake" collapse. Thomas's spherically 
symmetric code allowed for a multi-phase fluid but no collisionless dark matter, while the code 
we develop here evolves a single fluid component and a single collisionless component. Shapiro 
& Struck-Marcell included a collisionless component, and they examined collapses with planar 
rather than spherical symmetry. Our treatment of radiative cooling is somewhat different from 
that adopted by these authors; in particular, we can include the influence of a photo-ionizing 
background on the abundances of ionic species. However, the largest differences are not in the 
codes but in the choice of problem, and the consequent choice of initial conditions. 

The geometry in our calculations is idealized, and one must therefore take care to keep 
their limitations in mind. Nonetheless, they can provide a valuable complement to their more 
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elaborate, 3- dimensional cousins because of their high resolution, their speed, and their relative 
simplicity. Three-dimensional hydrodynamic simulations of galaxy formation suffer from limited 
spatial resolution and mass resolution, making it difficult to separate genuine physical results from 
numerical artifacts. One-dimensional collapse calculations can achieve much higher resolution, 
computing gas dynamic processes with much higher accuracy. Because they run fast, it is possible 
to undertake a much more comprehensive exploration of the parameter space, varying both the 
values of cosmological parameters and the assumptions about the gas microphysics. The results 
are much easier to visualize and interpret than those of three-dimensional simulations. Thus, these 
simplified, high-resolution calculations can provide a useful numerical check on three-dimensional 
simulations and, equally important, provide physical insight into their results. They can also 
check and improve upon the simpler analytic models that serve as inputs to Press-Schechter type 
calculations. 

This paper serves two purposes. First, it describes the code itself, and tests its ability to 
reproduce known analytic results such as the self-similar, spherical infall solution of Bertschinger 
(1985). Second, it applies the code to one of the basic questions of galaxy formation: what 
causes the abrupt cutoff at the upper end of the galaxy luminosity function? It has long been 
recognized that gravitational effects alone cannot explain this cutoff, because the largest virialized 
objects - rich galaxy clusters - have much higher masses than the largest galaxies (White & 
Rees 1978). The "lore", deriving largely from Binney (1977), Silk (1977), and especially RO, is 
that the upper cutoff is determined mainly by atomic physics, specifically by the requirement 
that the gas within a density perturbation be able to cool and collapse within a Hubble time. 
RO even include a "numerological digression" in which they relate the characteristic masses 
and sizes of galaxies directly to fundamental gravitational and atomic constants, independent of 
cosmological parameters. We will examine the underpinnings of this argument by studying the 
dynamics of gas in collapsing systems of various masses, focusing especially on the cooling in 
high-mass perturbations. In a later paper, we will examine the suggestion of Efstathiou (1992) 
that photo-ionization by the UV background may strongly affect the formation of low-mass 
galaxies. 

This paper is organized as follows. In §2 we present our general numerical model, focusing on 
the treatment of hydrodynamics and cooling in the gas component. In §3, we show the results of 
several test calculations. In §4, we present results for collapses of spherical density pertubations, 
with and without a collisionless component. In §5, we discuss the implications of these results for 
the galaxy luminosity function, comparing our analysis to those of RO and White & Frenk (1991). 

2. Numerical Model 

The simplest model for galaxy formation consists in the evolution of a uniform, pressureless, 
spherical density enhancement in a Friedmann Universe. The expansion of such a region lags 
behind the Hubble flow, until it stops at a turnaround time t ta and radius r ta , and recollapses. 
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It then undergoes violent relaxation and virializes after another ~ 1 — 2t ta - The value of t ta 
depends on the amplitude of the density contrast 6p/p at the recombination epoch. The large-scale 
distribution of galaxies is consistent with smooth hierarchical clustering resulting from a purely 
gravitational process, with no preferred length-scale. However, dissipation must have played a role 
in the formation of galaxies themselves. Within the spherical model, one expects gas to shock and 
heat as it collapses, and pressure forces and radiative cooling will strongly affect the post-shock 
evolution. 

To better understand these important effects, we have developed a simple but highly accurate 
one- dimensional numerical code to model the collapse of individual spherical perturbations. 
The code treats a mixture of gas particles, evolved through Lagrangian hydrodynamics, and 
collisionless particles (cold dark matter), each represented by concentric shells. The Lagrangian 
description is preferable to an Eulerian description because it follows the mass and fluid elements 
themselves, thereby maintaining the mass resolution throughout the calculation. As we will show 
in §4, the density and temperature are extremely non-uniform during the collapse, leading to 
widely varying timescales, both in space and in time. 



2.1. Equations 



Since we study isolated, collapsing, density perturbations, we will work with physical 
coordinates (rather than coordinates that are comoving with the expanding universe). 

The gaseous component is described by the fluid equations for a perfect gas. These are the 
continuity equation, 

^■p fl V-v g = > (1) 



the momentum equation, 
the energy equation, 
and the equation of state, 



dt 

g, (2) 



c?v g Vp 



dt Pg 

V - A 

(3) 



du p dp g T — A 



dt pj dt p g 



P = (7 " l)/> a «- (4) 



In these equations, p g , v g , p, and u are the baryonic mass density, velocity, pressure, and internal 
energy per unit mass, T is the external heating rate (e.g. from photo-ionization), A is the radiative 
cooling rate, 7 is the adiabatic index, and g = — V$, where $ is the total (gas and dark matter) 
gravitational potential. 

In spherical symmetry, equations (1) and (2) can be rewritten in the Lagrangian form as 

dm g = \-Kr 2 pgdr g , (5) 
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which replaces the continuity equation, and 

where r g is the radius of the gas shell and M(r) is the total (baryonic and dark matter) mass 
inside radius r. Here and throughout this paper we have set G = 1. 

The collisionless component is simply described by the equation of motion 

where is the dark matter velocity. In spherical symmetry, this can be rewritten as 

dv d M(r d ) 



dt rj ' 



(8) 



where rj is the radius of the collisionless mass shell. 



2.2. Radiative cooling 



We compute radiative cooling for a gas of primordial composition, 76% hydrogen and 24% 
helium by mass. The full set of equations that we use to obtain abundances and cooling rates 
is listed in Katz, Weinberg & Hernquist (in preparation), so here we restrict ourselves to a brief 
summary of the physics. Our original source for most of these formulae is Black (1981), and we 
adopt the high-temperature corrections of Cen (1992). 

We compute the abundances of ionic species as a function of density and temperature by 
assuming that the gas is in ionization equilibrium with a spatially uniform background of UV 
radiation. In other words, we choose the abundances so that the rate at which each species is 
depopulated by photo-ionization, collisional ionization, or recombination to a less ionized state 
is equal to the rate at which it is populated by recombination from a more ionized state or by 
photo-ionization or collisional ionization of a less ionized state. The intensity and spectrum of the 
UV background are specified as a function of time by the user, based on theoretical models or 
observational constraints. Given the ionic abundances and the density, we compute the cooling 
rates due to collisional excitation, collisional ionization, recombination, and Bremsstrahlung, and 
we compute the heating rate from photo-ionization. For the physical problems that we study here, 
the timescales for reaching ionization equilibrium are much shorter than the other timescales of 
interest, so our equilibrium assumption should be an excellent approximation. 

For the simulations in §§4.2 and 4.3 of this paper, we set the UV radiation background to 
zero. In this case the relative abundances are determined by collisional equilibrium alone, and they 
depend only on temperature. Since all of the cooling processes that we consider involve two-body 
interactions, we can describe the cooling rate A by a single function of temperature, up to a factor 
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Fig. 1. — Cooling rates as a function of temperature for gas of primordial composition (/# = 
n H m p/Pg = 0.76, fue = 0.24). The thick solid line is the total cooling rate. The dotted lines 
are the cooling rates from collisional ionization of H, He, and He + ; the long-dashed lines are 
the cooling rates from recombination to H, He, and He + ; the short-dashed lines are the cooling 
rates from collisional excitation of H, He, and He + ; the thin solid line is the cooling rate from 
Bremsstrahlung. Abundances of ionic species are computed assuming collisional equilibrium. 
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of Pg. The thick solid line in Figure 1 shows the total cooling rate A/n^, where nu = 0.76p g /m p is 
the number density of hydrogen nuclei. Short-dashed lines show the contributions from collisional 
excitation of H and He + , which dominate cooling in the range 10 4 K < T < 10 5 ' 5 K. Dotted lines 
show the contributions from collisional ionization. Long-dashed lines show the contribution from 
recombination. The thin solid line shows the Bremmstrahlung contribution, which dominates at 
T > 10 5 ' 5 K. Below 10 4 K, all of the gas is neutral, and there are no collisions energetic enough to 
cause electronic excitations, so none of the processes that we consider can produce any significant 
cooling in this regime. At very high densities the formation of hydrogen molecules can cool 
primordial gas to the temperatures required for fragmentation into stars. We will not attempt to 
resolve stellar mass scales in our cosmological studies, so we will just consider gas that cools to 
10 4 K to be "cold" and leave it at that. 

The principal scientific concern of this paper is the cutoff at the high-luminosity end of the 
galaxy luminosity function. We are therefore interested primarily in the behavior of high-mass 
perturbations. At the virial temperatures associated with these perturbations, typically 10 6 K 
or greater, the gas is fully ionized by collisional processes, so our neglect of the photo-ionizing 
background should make no difference to our conclusions. Photo-ionization can affect the behavior 
of lower mass perturbations because it eliminates the neutral hydrogen and singly ionized helium 
that dominate cooling at low temperatures, and because the residual energy of the photo-electrons 
heats low-density gas to T ~ 10 4 K. Our next paper will focus on the influence of photo-ionization 
on the collapse and cooling of low-mass perturbations. 

Compton scattering of microwave background photons by electrons can be an important 
source of cooling at redshifts z <; 10. However, in the collapse calculations in this paper the gas 
remains neutral until much lower redshifts, when it collapses and shocks, so we can safely ignore 
Compton cooling. It would be straightforward to add this effect to our code, and we would need 
to do so in order to study collapses at higher redshifts. 

2.3. Numerical scheme 

We use the standard, second-order accurate, Lagrangian finite-difference scheme (Bowers 
& Wilson 1991). In this scheme, the velocity is zone-edge-centered, while the pressure and 
internal energy are zone-centered. To obtain time-centering, the velocities are evaluated at 
half-timesteps. We give the hydrodynamical finite-difference equations in the order in which 
they must be evaluated. In the following, the subscripts denote the position of the shell, and 
the superscripts denote the time. We first present the equations for the gas component, and for 
clarity of presentation we drop the subscripts g. Note that in the following equations, p™ is the gas 
density, and m™ is the total mass interior to the shell at position i at time n, including both gas 
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and dark matter. First we must advance the velocities to f n+1 / 2 according to 

dt n . (9) 



n+l/2 n-1/2 
V- = V- 



a I n^2 P i+l/2 Pi-1/2 , m T - 



drrii (rf) 2 
Then we can advance the positions and evaluate the densities, 

r^ 1 =r? + ^ +1/2 d* n+1 / 2 , (10) 



and 



In these equations, 



= d ^+l/2 (n) 



and 



dt n = ^(dt 71 - 1 / 2 + dt n+1 l 2 ), (12) 

dm, = ^(c?m l _ 1/2 + dm l+1 / 2 ). (13) 

Equations (9)-(ll) are second-order accurate in space and time. We can now advance the energy 
equation: 

„n+l n n [ l _ l _\ , ^ ~ A )»+l/2 , n+l/2 /^x 

U i+l/2 - U i+l/2 Pi+1/2 n+1 „n + n+1 <" • \ L *J 

\^+l/2 n+l/2 J Pg,i+l/2 

This last equation is only first-order accurate in time. To make it second-order accurate, we would 

have to replace 2>™ +1 / 2 an( ^ — -^)i+i/2 Pi+1/2 an( ^ — ^-Yi+i/2' This would require two 
evaluations of the cooling function per timestep, and since a lot of the computational time is in 
practice spent evaluating the cooling functions, this would be quite expensive. However, as we will 
show in §3, the energy is extremely well conserved in our present scheme, with equation (14), and 
there is therefore no reason to require more accuracy in the energy equation. 

Shocks are treated with the usual artificial viscosity technique (Richtmyer & Morton 1967). 
The pressure in the momentum and energy equations is replaced by P = p + q, where 



"+i - _ c - \ v n+1/2 - v n+1 / 2 \( v n+1 / 2 - v n+1/2 ) (15) 

^+1,2- C \H P ^ I2 ) H II P1+ J V ^ } (15) 



if Vi+\^ 2 — v™ +1 ^ 2 < and q = otherwise. We use c q = 4, which spreads shock fronts over 4-5 
shells. 

The second-order accurate, finite-difference equations for the collisionless equation of motion 
are much simpler than those for the gas, 



n+l/2 _ n-1/2 _ ^ , lg x 

V d,i ~ V d,z (r£.) 2 ' ^ ' 



and 



^,r=^+^T /2 ^ +i/2 - (17) 
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In the absence of shell-crossing, the mass inside a given shell would be constant in time. 
However, the collisionless shells are allowed to cross each other and to cross gas shells, so the 
masses m™ are functions of time and must be computed at each time step. 



2.4. Timesteps and central boundary conditions 



For each shell in the calculation there are three potentially important timescales, namely the 
dynamical, Courant, and cooling timescales. In addition to respecting these timescales, we must 
ensure that the fluid shells do not cross. We therefore take a timestep 



dt = min {dtjyn, die our-, dt coo i-, dt ve i} , 



where 



dt, 



dyn 



mm 



dt 



Cour 



mm 



dt, 



cool 




and 



m i c v } , 

l I Vi- v^x J 



mm 



(18) 

(19) 

(20) 
(21) 

(22) 



where c^, cc, c c , and c„ are safety constants. We use = 0.01, cq = 0.2, c c = 0.1, and c„ = 0.05. 
The timesteps vary widely both in space and in time during a calculation. It is therefore essential 
to compute appropriate timesteps for all shells in the calculation at a given time, and to use the 
smallest of these to advance the system. We could increase computational efficiency at the price 
of additional complication by allowing different shells to have different timesteps, but thus far we 
have not found it necessary to adopt this procedure. 

A Lagrangian code achieves very high spatial resolution near the center, and it can thereby 
demand extremely short timesteps. For collisionless particles th e times tep is determined by 
dtdyn oc \/r 3 /m(r). The first shell, with mass dm, has dtd yn oc ^r\/dm oc r^ 2 , so the timestep 
goes to zero as the shell approaches r = 0. Since spherical symmetry is an idealization in any 
case, there is no need to bring the calculation to a grinding halt in order to integrate the very 
central region at high accuracy. We therefore solve the timestep problem by treating the center as 
a hard reflecting sphere of radius r c (Spitzer & Hart 1971; Gott 1975), which prevents timesteps 
from becoming arbitrarily small. Energy conservation degrades as r c is decreased (because shells 
rebound at higher velocities), but it is not difficult to find a value of r c that (a) yields excellent 
energy conservation and (b) is much smaller than other characteristic radii in the problem, so that 
the departure from an idealized spherical collapse is minimal. 
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A similar problem arises for the gaseous component when the gas elements cool very rapidly. 
In regions where the cooling time is much smaller than the dynamical time, the fluid shells 
collapse at the free-fall rate, the central density increases very rapidly, and the cooling and 
dynamical timesteps become exceedingly small. We solve the problem by "freezing" shells that 
have t coo i <C tdyn- More specifically, when a shell has t coo i < Cftd yn , we move it to a radius 
r = r m i n , assign it zero velocity and a temperature of 10 4 K, and subsequently ignore it, except in 
the calculation of the gravitational force. We require (r — r m i n ) < 0.1r m i n , so that no shell moves 
a large fraction of its radius in the freezing process; if the condition t coo i < Cftd yn is reached at 
a large radius, we continue to evolve the shell until it reaches r < \.\r m { n . Shells in this regime 
collapse nearly isothermally at the free-fall rate. We have run tests with Cf = 10~ 2 and Cf = 10~ 3 ; 
these values yield virtually identical results. We adopt r m i n = r c , the radius of the reflecting 
sphere used for the collisionless component. We first tried moving frozen shells to the origin, but 
this practice leads to unstable numerical results because the shells have a wide range of values of 
tdyn at the time they are frozen. Adopting r m i n = r c suppresses this instability, though in one 
physical regime (where a large fraction but not 100% of a perturbation collapses), the amount 
of cooled mass can vary by ~ 10% from one run to another, depending on the number of shells 
used, because of residual sensitivity to the central boundary condition. Such variations are small 
compared to the approximation made in treating galaxy formation as a spherically symmetric 
process. 

3. Test Calculations 

We have performed a variety of tests of the code on problems with known analytic solutions. 
We describe results from several of these tests in this section. 

3.1. Polytropes 

In the absence of radiative cooling, the total energy of the system must be conserved, and 
in the absence of shocks, the entropy must be conserved. To check that the code satisfies these 
basic requirements, we simulate an equilibrium polytrope of index n p = 1.5 (7 = 5/3). We use 
N g = 500 shells, equally spaced in radius r, and units such that G = M = R = 1. We relax the 
polytrope to an equilibrium by adding a dissipative term in the momentum equation of the form 
—Vg/Uelax- We then remove this term and allow the polytrope to evolve dynamically. The result 
is an exceedingly accurate equilibrium: positions of the fluid shells fluctuate with amplitudes 
Ar/r smaller than 10~ 7 , the kinetic energy remains smaller than 10~ 14 , and the potential energy, 
thermal energy, and entropy change by less than 10~ 8 over ~ 100 dynamical times. 

For a more rigorous trial, we set up initial conditions corresponding to the first two 
normal pulsation modes of the polytrope. To obtain the appropriate initial configuration, we 
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Fig. 2. — (a) Eigenfunctions for the first two pulsating modes of a polytrope with adiabatic exponent 
j = 5/3. To set up a pulsating polytrope we perturb the shell radii of the equilibrium model 
according to r — > r + 8. (b) Fluctuations in the total energy, the gravitational potential energy, the 
internal energy, and the kinetic energy, from a simulation with N g = 500 shells, during 20 periods 
of the second normal mode. The bottom panel shows entropy production during the oscillations, 
(c) Fluid shell trajectories and entropy production when the second normal mode has an amplitude 
of 15% of the star radius. 
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simultaneously integrate the Lane-Emden equation and the eigenvalue equation governing the 
radial oscillations of a spherical star (Shapiro & Teukolsky 1983). We choose an initial mode 
amplitude equal to 1% of the star radius. The eigenfunctions are shown in Figure 2a. The 
eigenfrequencies corresponding to these eigenfunctions are oj\ = 1.645, and u>2 = 3.547. We let the 
pulsating polytrope evolve dynamically, using N g = 500 shells in each case. The total energy and 
entropy vary by less than 5 X 10~ 5 over about 20 periods of oscillation. The oscillations of the 
gravitational, thermal, and kinetic energies are shown in Figure 2b. The oscillation periods are 
T\ = 3.82 and T2 = 1.77 for the first and second modes, respectively. These are identical to the 
periods 2tt/u> obtained from the eigenfrequencies of these modes. 

From Figure 2b we see that the entropy grows slowly during evolution. This entropy increase 
arises because of the finite amplitude of the modes. If we increase the mode amplitude further, 
the oscillations produce shocks, which generate entropy. Figure 2c shows this entropy generation 
when the second normal mode has an amplitude of 15% of the star radius. The shock appears in 
the outer shells, where the amplitude of the perturbation is the largest. 



3.2. Self-similar hydrodynamic collapse 



To test the code in a context closer to its cosmological purpose, we check that it reproduces 
the similarity solution for shocked accretion of a collisional, non-radiative gas, as described by 
Bertschinger (1985). This solution describes the asymptotic (late-time) behavior of a spherically 
symmetric collapse about an initial seed perturbation in an Einstein-de Sitter (f2 = 1) universe. 
In the absence of cooling, neither the gas physics nor the cosmology defines a preferred scale, so 
once the collapsed mass exceeds that in the initial perturbation, the system "forgets" the details 
of its initial state, and its evolution becomes self-similar in time. We evolve a perturbation with a 
Gaussian initial density profile, 

6 l (r) = 6 l {r = 0)e- r2 l R *, (23) 
and an unperturbed, Hubble-flow velocity profile, 

Vi(r) = H.r, Hi = 2/(3*0- ( 24 ) 
Here t{ is the initial time, H{ is the Hubble constant, and 

fc(r)=[*(r)-a]/ft (25) 

is the initial density contrast, with 

■p l = Pm = {&-Ktl)-\ (26) 
the critical density at time t{. In the linear regime, the density contrast grows as 
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where 

— . x f n r 6(r')4Trr' 2 dr' , x 

6(r) = J ° V 3 (28) 



3 -ir 



is the averaged overdensity interior to the shell at position r, and 

r = t/U (29) 

is the time in units of the initial time t{. The second term in equation (27), and the 3/5 factor 
multiplying the first term, appear because our Hubble-flow initial conditions contain a mixture of 
growing and decaying modes. 

As the density contrast inside a shell grows, the extra gravitational deceleration drags it 
further behind the Hubble flow, until it finally turns around and recollapses. Roughly halfway 
back to the center, it hits a shock, which sharply raises its density, temperature, and entropy, and 
brings it nearly to rest. Thereafter, the shell falls very slowly towards r = 0. Figure 3 shows 
the trajectory of a typical shell in our simulation; the radius and time are scaled to the shell's 
turnaround radius r' ta and turnaround time t ta - With this scaling, the trajectories of all shells 
that lie well outside the initial perturbation are identical, and the trajectory shown in Figure 3 is 
indistinguishable from that in figure 4 of Bertschinger (1985). 

Figure 4 shows the velocity, density, pressure, and mass profiles from a simulation with 
<^(r = 0) = 0.2, Ri = 1, a 7 = 5/3 equation of state, and N g = 1000 shells. Solid lines show results 
at r = 1000, 2000, 3000, 4000, and 5000, from bottom to top. The dashed line, often obscured by 
the solid lines, shows Bertschinger's (1985) similarity solution. The dimensionless radius, velocity, 
mass, density, and pressure are defined by 



A = r/r ta , 


(30) 


V{\) = v{r,t){r ta lt)-\ 


(31) 


M(\) = m(r,t)(^Tr PH r^ , 


(32) 


D(X) = p(r,t)/p H , 


(33) 


P(\)=p(r,t)(t/r ta ) 2 pJ I 1 . 


(34) 



and 



Here r ta is the radius of the shell that is currently turning around (and is thus distinct from 
the radius r' ta used in Figure 3). We see in Figure 4 that the profiles for the dimensionless fluid 
parameters tend asymptotically towards the similarity solution, demonstrating the ability of our 
code to reproduce known analytic results for this problem, one that is directly relevant to the 
cosmological applications that we will consider. The total energy is conserved to better than 0.2% 
over the entire run. 



- 15 - 




Fig. 3. — Trajectory of a fluid shell in the self- similar, shocked accretion of a 7 = 5/3 collisional 
gas. The radius and time are scaled to the shell's turnaround radius r' ta and turnaround time t ta - 
With this scaling, the trajectories are identical for all shells. 
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Fig. 4. — Convergence of the numerical results toward the similarity solution for shocked accretion 
of a 7 = 5/3 collisional gas. The non-dimensional velocity V, density D, pressure P, and mass M 
are shown as a function of the scaled radius A. The units are given by eqns. (30)-(34). Profiles 
from a simulation with N g = 1000 shells are plotted at r = 1000, 2000, 3000, 4000, and 5000, from 
bottom to top, where r is defined in eqn. (29). The dashed line shows the similarity solution from 
Bertschinger (1985). 
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3.3. Pressureless collapse onto a black hole 

In the case of cold accretion onto a black hole, there is no shell-crossing, M(r) is constant 
in time, and the equation of motion r = — M(r)/r 2 can be integrated analytically for each shell. 
Integrating once we get v 2 /2 = M/r — C, where C = M/ri — vf/2, and T{ and V{ are the position 
and velocity of the shell at the initial time t{. This equation can be integrated once more to 
give the radii as function of time, yielding the well-known cycloid solution, usually written in 
parametric form (see, e.g., Padmanabhan 1993, §8). The mass profile and the velocity profile are 
therefore known analytically. 

For this test, we start with a top-hat initial density profile, 

(1 + 6, iir<R l 
" = '"*il ifr>*, (35) 

with Si = 0.3 and Ri = 1. We evolve the collisionless system forward in time and compare the 
numerical results for the mass and velocity profiles to the exact analytical values. In Figure 5a 
and 5b, solid lines show the numerical velocity and mass profiles, obtained for Nj = 10,000 shells. 
Points show the exact results at selected values of A. Dashed lines represent the unperturbed 
Hubble flow. Numerical and analytic results agree to better than 0.5%, as shown in Figures 5c 
and 5d. 

At the time shown in Figure 5, the collapsed mass significantly exceeds the mass in the 
initial top-hat perturbation, so the system has reached the regime of self-similar evolution. 
Our Figures 5a and 5b are, therefore, identical to the similarity solution plotted in figure 1 of 
Bertschinger (1985). 

3.4. Coupled hydrodynamic and collisionless systems 

One important new numerical effect enters when we model mixed collapses: as collisionless 
shells cross fluid shells they introduce discrete fluctuations in the gravitational forces acting on 
the fluid, and this spurious agitation leads to low-level shocks and associated entropy production. 
The discreteness effect depends primarily on the mass ratio between collisionless and fluid shells. 
From a variety of tests, we find that this effect becomes negligible when the mass ratio is unity or 
smaller. In general, therefore, we require Nj/Ng > Qj/Qg. 

4. Collapses with radiative cooling 

We now turn to the main scientific issue of this paper: the relation between the physics of gas 
cooling and the rather sharp upper cutoff in the distribution of galaxy luminosities. We address 
this point by modeling spherical collapses with radiative cooling appropriate to a gas of primordial 
composition in collisional equilibrium (see §2.2 and Figure 1). 
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Fig. 5. — Numerical results (solid lines) for pressureless collapse onto a black hole, from a simulation 
with Nd = 10,000 collisionless shells. Upper panels show the non-dimensional velocity V and mass 
M as a function of the scaled radius A. The units are given by eqs. (30)-(33). Points show exact 
analytic results, and dashed lines show the result for unperturbed Hubble flow. Lower panels show 
the error in the numerical results. The velocity error is scaled to a characteristic velocity of the 
problem, vh(^ = !)• 
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4.1. Initial conditions 



As initial conditions, we adopt the average density profile around a peak in a Gaussian 
random density field. Equation (7.10) of Bardeen et al. (1986) gives the orientation- averaged, 
mean density profile around a peak of height vo{ and curvature x, 

F(r) v fli x/P fa2l , 

— = J^W)^ + ~ } " o=W)W * + ~ } ' (36) 

where /3 is a constant that depends on the slope of the power spectrum, ipir) = £(?")/£(0) is the 
normalized correlation function, and Oi is the rms density fluctuation. For a density field with a 
power-law power spectrum of index n smoothed by convolution with a Gaussian filter of radius 
Rf, the power spectrum is P(k) = Ak n e f , and the normalized correlation function, its Fourier 
transform, is 

*> = ^!ii}>' < 37 > 

where </> is the degenerate hypergeometric function (defined in Gradshteyn & Ryzhik 1980, 
§9.210.1). In this paper we adopt n = —2, roughly the slope of the cold dark matter power 
spectrum on galactic scales. For this choice, </> is related to the gamma function through 
7(a, x) = (x a I a)</>(a, 1 + a; —x). The coefficient /3 = \Z(n + 3)/(ti + 5) = 1/3, and 

3(1 - (5 2 ) + (1.216 - 0.9/3 4 )e-^/W 2 ) 2 
<x>-pv+ [3(1 _ ^ + Q A5 + (/3z//2)2]1/2 + p vj2 ■ (38) 

Figure 6 shows the normalized profile, F(r)/ai, and the spherically averaged interior overdensity, 
F(r)/(Ti, for a loi peak. It is this profile that we adopt as the initial density distribution for our 
simulations. We do not expect our qualitative results to be sensitive to the details of this choice. 
The perturbed initial velocity profile is that implied by the growing mode solution of linear theory, 
Vi = Hiri(l - 6i/3). 

Given the profile shape of Figure 6, our initial conditions have two free parameters: the filter 
radius Rf and the initial overdensity of the peak at r = 0, 6 p . Physically, it is more convenient 
to describe the perturbations in terms of a mass scale and a redshift of collapse. We define the 
filter mass Mf to be the mass contained within a sphere of radius 2Rf, and we define the collapse 
redshift z c to be the redshift at which the r = 2Rf shell would collapse to r = in the absence 
of pressure. For our adopted profile, 6(2Rf)/6 p = 0.6805, so the mass Mf is related to the filter 
radius Rf by 

M f =f(2Rf) 3 p Htl [l + 8 l (2R f )] (39) 
= ^ f (2 J ff / ) 3 (l + 0.6805<y. (40) 

i 

Specifying the initial time and the collapse redshift z c determines the value of 6 P , 

6 P = 6i(2R f )/0.6805 = 1.69 (^-^-\ /0.6805. (41) 

V 1 + Zi J 
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Fig. 6. — Mean profile of a 2<r peak in a Gaussian density field with power spectrum P(k) = 
Ak~ 2 e~ k R f . The solid line shows the density contrast; the dashed line shows the mean interior 
density contrast, as defined in eqn. (28). We adopt this profile as the initial density distribution 
for our collapse simulations. 
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Equation (41) is specific to an = 1 universe, but the generalization to an open universe is 
straightforward. 

By varying the mass Mf collapsing at a fixed redshift z c , we effectively vary the filter radius 
Rf. By varying the redshift at which a given mass collapses, we vary the density at which the 
collapse occurs. In practice, we want to ensure that simulations start in the linear regime, so we 
fix the initial peak overdensity at 8 V = 0.2. Specifying z c determines the initial redshift through 
equation (41), and specifying Mf determines Rf through equation (40). The initial redshift is 
Zi ~ 36 for z c = 2 and Z{ ~ 24 for z c = 1. 

4.2. Pure fluid collapses 

In the absence of cooling, physical processes do not introduce any preferred length or time 
scales. In appropriate units, therefore, the results of a collapse are the same for all values of Mf 
and z c : radii scale as r oc M^ 3 (l + -z c ) _1 , and times scale as t oc (1 + z c )~ 3 / 2 . Figure 7 shows the 
r vs. t trajectories of several fluid shells in a collapse calculation without cooling, using the initial 
density profile described in §4.1 and N g = 300 fluid shells. Shell radii are scaled to the initial 
filter radius Rf, and times are scaled to the collapse time t c of the shell with initial radius 2Rf. 
In these units, the fluid trajectories are identical for any values of Mf and z c . Each shell expands 
initially, then turns around and recollapses until it is halted by a shock, usually when it has fallen 
back to about half of its maximum radius. The shock forms at the center when the first shells 
collapse, and it propagates outward in both Eulerian and Lagrangian coordinates. When a shell 
hits the shock, its fluid density increases (by a factor of four in the limit of a strong shock), and 
nearly all of its kinetic energy is converted to thermal energy. At r s = r ta /2, the pre-shock infall 
velocity is v = (GM/^) 1 / 2 , making the post-shock temperature T ~ (GM /r s )(p,m p /3k), roughly 
the virial temperature implied by the interior mass M and the shock radius r s . After the shock, 
the shell is almost supported by the pressure of the hot gas beneath it, but as more gas piles on 
top it is compressed very slowly towards r = 0. The trajectories of the outer shells in Figure 7 
are the same as the shell trajectory for the similarity solution (Figure 3), but their shape appears 
somewhat different because they are plotted in different units. 

Radiative cooling introduces a new timescale into the collapse problem. Prior to shocking, 
pressure is unimportant, so a shell follows the same trajectory as before. However, the behavior 
of the shell after it hits the shock depends critically on the post-shock density and temperature, 
specifically on the ratio of the cooling time, t coo i ~ u p/A-, to the dynamical time, td yn ~ (G/o) -1 / 2 , 
where p is the mean mass density interior to the shell. While the local density p and the 
mean interior density p are not identical, they are roughly proportional to each other. At fixed 
temperature, the cooling rate A oc p 2 (assuming collisional equilibrium), so t coo i/td yn ~ p~ X ^ 2 ■ 
A shell that collapses earlier, when the density is higher, cools more efficiently. The post-shock 
temperature determines the initial location of the shell on the cooling curve (Figure 1), and 
because the cooling curve is a complicated function of temperature, the influence of the post-shock 



- 22 - 




12 3 



t/t c 

Fig. 7. — Shell trajectories for pure fluid (7 = 5/3) collapses, in the absence of cooling. The time 
and radii are scaled with respect to the collapse time t c and the filter radius Rf, defined in §4.1. 
With this scaling, the trajectories are independent of Mf and z c . 
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temperature on shell behavior is itself rather complicated. In particular, it is important to recall 
that the radiative cooling cuts off very sharply at T ~ 10 4 K. 

We have performed a series of collapse calculations in which we vary the value of the filter 
mass Mf while keeping the collapse redshift fixed at z c = 2. Each calculation uses N g = 1000 
shells. The minimum radius (see §2.4) is set to 0.99 times the initial radius of the innermost shell, 
making r c = OARf. We performed tests where we used N g = 5000 shells and r c = 0.2Rf. The 
results were essentially identical. 

Figures 8-10 illustrate four quite different histories of individual shells from these calculations, 
with Figure 8 showing the r vs. t trajectories, Figure 9 the trajectories in the n — T plane, and 
Figure 10 the time evolution of the temperature and the timescale ratio t coo i/td yn . In Figures 8 and 
10, to is the age of the universe at z = 0. The collapse redshift z c = 2 corresponds to t = 0.192io- 
In Figure 9, nn is the mean hydrogen number density inside the shell radius, nn = fHf>/ m p, with 
f H = 0.76 the mass fractio n of hydrogen. The dotted lines in this Figure show contours on which 
the ratio t coo i/td yn is constant, demarcating regions where cooling is rapid or slow relative to the 
dynamical time. The structure of these contours is closely related to the structure of the cooling 
curve itself. Dashed lines indicate contours in the n — T plane along which the cooling time or the 
dynamical time is equal to the age of the universe. Diagonal solid lines are lines of constant Jeans 
mass, Mj = {itk B I Gm^fl 2 ^! 2 p- 1 ! 2 . 

Because Figure 9 plots the quantity Hh corresponding to the mean interior overdensity, which 
is the relevant parameter for the dynamical time, the computation of t coo i, which depends on 
the local density, is somewhat ambiguous. In the similarity solution the ratio of local density in 
a recently shocked shell to the mean interior density is about a factor of four. With this result 
in mind, we have computed t coo i(nH,T) in Figure 9 on the assumption nn = 4tTh- While this 
approximation is better than nn = nn, one should still take the contours in Figure 9 as indicative 
rather than precise boundaries. The cooling times in Figure 10 are computed from the shell's local 
density, so the t coo i/td yn ratios plotted there are more reliable. 

Our four illustrative fluid shells are labeled A-D in these Figures. In Figure 8 the dot-dash 
line shows the trajectory of a shell with no cooling for comparison. Shell A has an interior mass 
of 6.3 X 10 7 M e = 2.9Mf. Its post-shock temperature is smaller than 10 4 K, so after the shock the 
shell is unable to cool, and it remains nearly pressure supported. However, the continuing infall 
of gas from larger radii compresses the shell, increasing its density and temperature adiabatically. 
The ratio of cooling time to dynamical time is large during this phase of evolution (Figure 10). 
Eventually, compression pushes the shell temperature above 10 4 K; the cooling time drops rapidly, 
and the shell collapses to r = at the free-fall rate. This phase of the collapse is effectively 
isothermal, since any energy gained during compression is immediately radiated away. The points 
on trajectory A in Figure 9 are spaced at equal time intervals At = 0.1 to. One can see from their 
locations that the initial post-shock evolution is very slow but the final collapse very rapid. All 
points above lognn = — 1 in Figure 9 belong to trajectory B; the next point for shell A lies off the 
top of the plot. 
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Fig. 8. — Trajectories of fluid shells in pure fluid collapses with radiative cooling, for collapse 
redshift z c = 2 and four different values of Mf. Dimensionless radii r/Rf are plotted as functions 
of the dimensionless time t/to, where to is the age of the universe at z = 0. The masses interior to 
the shells are6.3xlO 7 M = 2.9M f (shell A), 3.2 x 1O 8 M = 3.2M f (shell B), 3.2 x 1O 14 M = 3.2M f 
(shell C), and 3.2 X 1O 16 M = 3.2M f (shell D). The dot-dashed line shows the shell trajectory from 
a calculation without cooling. 
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Fig. 9. — Post-shock trajectories in the log(T)-log(TC/j) plane of the four fluid shells shown in 
Figure 8, labeled A-D as before. The trajectory of A is shown by a dashed line, and the trajectories 
of other shells by solid lines. Points are spaced at intervals At = O.lio along trajectory A, 
At = 0.0012£o along trajectory B, At = 0.025io along trajectory C, and At = O.lio along trajectory 
D. Dashed lines show the contours t coo i = to and td yn = to. Dotted lines represent contours of 
constant t coo i/td yn . Cooling times are not precise because they depend on local density nn instead 
of mean interior density nn', we compute them assuming nn = ^nn- Diagonal solid lines are 
contours of constant Jeans mass Mj. 



- 26 - 




Fig. 10. — Time evolution of the temperature (upper panel) and the ratio t coo i/td yn (lower panel) 
for the four fluid shells shown in Figures 8 and 9. 
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Fluid shell B has an interior mass of 3.2 X 10 M = 3.2Mf. For this shell the post-shock 
temperature is just slightly higher than 10 4 K, in a regime where cooling is very rapid. The shell 
collapses to r = isothermally, at the free-fall rate, with no phase of pressure support. The 
post-shock evolution is extremely rapid; points along trajectory B in Figure 9 are evenly spaced in 
time, but the interval is only At = 0.0012 to. 

Shell C has an interior mass of 3.2 X 1O 14 M0 = 3.2M^, six orders of magnitude higher than 
that of shell B. Immediately after the shock, the cooling time exceeds the dynamical time by more 
than a factor of 10. Therefore, unlike shell B, shell C goes through a period of nearly adiabatic 
compression, with little radiative cooling. During this phase, the shell stays in quasi-static 
equilibrium close to the border of Jeans stability, so it follows a track of constant Jeans mass in 
the n — T plane (Figure 9). As the density increases, the ratio of cooling time to dynamical time 
decreases steadily. Once t coo i < td yn , the shell cools rapidly to 10 4 K and collapses isothermally at 
the free-fall rate. Points along trajectory C in Figure 9 are spaced at intervals At = 0.025 to. 

Shell D has an interior mass of 3.2 X 10 16 = 3.2Mf. The post-shock cooling time is again 
much longer than a dynamical time. In fact, the cooling time in this case is comparable to the 
age of the universe, so shell D evolves quasi- statically and never enters a rapid cooling phase. If 
the simulation were evolved further, it would eventually reach a density high enough for cooling 
to become important, and it would behave in a fashion more similar to shell C. Points along 
trajectory D in Figure 9 are spaced at intervals At = O.lto. 

If we want to understand features of the galaxy luminosity function, the quantity of most 
interest is the mass of gas that actually cools by redshift zero, since this is the gas that could 
potentially fragment into stars. The filled circles in Figure 11 show the ratio M c /M s as a function 
of log(M s ), where M c is the mass of gas that collapses, shocks, and cools by z = 0, and M s is the 
mass of gas that collapses and shocks by z = 0. For M s Ss 10 8 ' 5 Mq, the post-shock gas remains 
below 10 4 K, so it never cools, and M c = 0. Between M s = 10 8 - 5 Mq and M s = 10 n MQ, virtually 
all of the shocked gas cools, and M c /M s ~ 1. Above M s = IO^Mq, the outermost shocked shells 
remain pressure supported all the way to z = and do not cool. However, the inner shells of these 
perturbations collapse earlier, at higher density and lower virial temperature, and these shells are 
able to cool. Therefore, even though the ratio M c /M s decreases with increasing M s in this regime, 
the decrease is quite slow. In particular, M c /M s falls much less rapidly than M" 1 , so the actual 
amount of cooled gas increases with M s . We will return to this important point in §5. The results 
in Figure 11 are not sensitive to the assumed collapse redshift; we have carried out the same series 
of calculations for z c = 1, and the trend of M c /M s versus M s is very similar. 

4.3. Mixed collapses 

The basic physical effects described above carry over to the case of a mixed collapse involving 
gas and collisionless dark matter. Gas shells still collapse and shock, and their subsequent behavior 
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Fig. 11. — Ratio of the mass M c of gas that cools by z = to the mass M s of gas that shocks 
by z = 0, as a function of the virialized mass M s /Qb. Filled circles represent pure fluid collapses 
(fib = 1) with collapse redshift z c = 2. Open circles represent mixed collapses with = 0.1, 
Qd = 0.9, and z c = 2. 
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again depends on the cooling time and dynamical time at the post-shock density and temperature. 
However, in the presence of a mixture of dark matter and gas, the cooling time depends only 
on the gas parameters, t coo i ~ up g /A, while the dynamical time depends on the total mass, 
tdyn ~ Ptoti where p tot oc M(r)/r 3 . If we assume that the gas and dark matter density profiles are 
similar until the point where the gas cools, the ratio between p g and p to t is simply given by Qb/Q. 

— 1/2 

Therefore, we have t C ooll^dyn Ptot ^/^bj an d the ratio of the cooling time to the dynamical 
time is larger than that in the pure fluid case by a factor 0/0&. As a consequence, we expect the 
transition at high masses between the region where all of the gas cools and the region where only 
part of the gas cools to occur at a lower total mass. 

We have again performed a series of collapse calculations in which we vary the value of the 
filter mass Mf while keeping the collapse redshift fixed at z c = 2. In these calculations, we use 
Qd = 0.9, Qb = 0.1, N g = 500 fluid shells, and Nj = 10,000 dark matter shells. The open circles 
in Figure 11 show the results for the ratio of cooled gas to shocked gas, M c /M s , as a function of 
log (M s / f^). Once again M c is the cooled gas mass at z = and M s is the shocked gas mass at 
z = 0. The mass M s /Qb is, roughly, the total virialized mass at z = 0. As expected, the transition 
at high masses occurs at a lower threshold than in the case of a pure fluid collapse, but the 
general shape of the curve is similar, and the conclusion is the same, i.e the mass that cools is a 
monotonically increasing function of the total mass of the perturbation. 

5. Discussion 

Figure 11 shows three distinct regimes: at very low masses, M s <, 10 8 - 5 M Q , there is no cooling 
at all, at intermediate masses, 10 8 ' 5 Mq Sa M s Ss 10 11 Mq, all of the shocked gas cools, and at high 
masses, M s <; 10 11 Mq, only a fraction of the shocked gas cools. However, the transition from the 
intermediate-mass regime to the high-mass regime is not a sharp one. Five orders of magnitude 
above the transition mass, the fraction of shocked gas that is able to cool is still ~ 20% in the pure 
fluid case and ~ 5% in the mixed fluid/dark matter case. We see, therefore, that the requirement 
that gas be able to cool within a Hubble time cannot by itself explain the sharp upper cutoff in 
the luminous mass of observed galaxies. This point is further illustrated in Figure 12, where we 
show the mass of gas that cools by z = 0, M c , as a function of the mass of gas that has been 
shocked, M s . The cooled mass M c increases monotonically with M s , and the transition between 
the intermediate- and high-mass regimes is marked only by a modest change of slope in the M c 
vs. M s relation. For large masses the virial temperature is high, the gas is collisionally ionized, 
and the cooling is dominated by free-free transitions. Including a photo-ionizing background 
would not, therefore, alter our results at high masses. The low-mass behavior could be sensitive 
to assumptions about photo-ionization (Efstathiou 1992), a point that we will address elsewhere. 

RO suggest that cooling requirements can explain the turnover in the galaxy luminosity 
function, but the numerical results in Figure 12 indicate otherwise. Larger mass collapses can 
always produce larger mass galaxies, even if they do so with imperfect efficiency. In essence, the 
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Fig. 12. — The mass M c of gas that cools by z = as a function of the virialized mass M s /Qb. 
Filled and open circles represent pure fluid and mixed collapses, respectively. Solid lines show 
M c = M s , the result that would be obtained for complete cooling of shocked gas. The numerical 
results fall below these lines at low masses (where there is no cooling) and at high masses (where 
there is partial cooling), but there is no sharp transition at high masses that might correspond to 
the turnover in the galaxy luminosity function. The dashed line shows the result of applying the 
White & Frenk (1991) analysis to the case of pure fluid collapses. Filled triangles show, for the 
pure fluid case, the mass that cools within a single post-shock dynamical time. 
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difference between our result and RO's is the difference between a multi-zone and a single-zone 
calculation. RO associate a single characteristic density and a single characteristic temperature 
with the collapse of a given perturbation. They then ask whether a cloud of gas at that density 
and temperature can cool within a Hubble time, or within a dynamical time. However, typical 
collapses produce peaked density profiles rather than uniform profiles, so the inner, high-density 
regions can cool more efficiently than the outer regions. This point has been made in a somewhat 
different guise by White & Frenk (1991), in their semi-analytic models of the galaxy formation 
process. They assume that a collapse without cooling would produce a cloud with an r~ 2 density 
profile and a temperature equal to the halo virial temperature. They then compute the "cooling 
radius" — the radius out to which gas is dense enough to cool within a Hubble time — and from 
this they compute the cooled mass. This line of reasoning leads to good qualitative agreement 
with our numerical arguments, as shown by the dashed line in Figure 12, which represents the 
result of applying the White & Frenk analysis to the case of pure fluid collapses. The transition 
between complete cooling and incomplete cooling is gradual and subtle, unlike the turnover in the 
galaxy luminosity function. 

Our spherically symmetric calculation models the process of galaxy formation with a coherent 
collapse. However, the complex and untidy assembly of a proto-galaxy in a hierarchical scenario 
still generates a wide range of densities and corresponding cooling times, since material that 
collapses early reaches higher density and cools more efficiently. We therefore expect that the 
qualitative trend in Figure 12 would continue to hold in a more realistic calculation. Small-scale 
clustering in a hierarchical model will also heat gas prior to the proto-galaxy collapse, but the 
agreement between our numerical results and the White & Frenk analysis (which assumes that 
all gas is at the virial temperature) implies that it is the range of densities rather than the range 
of post-shock temperatures that explains the cooling in high-mass collapses. We have explicitly 
checked this point in the case of a pure fluid, Mf = 1O 15 M0 collapse, by running a series of 
simulations in which the gas is held at a finite temperature prior to shocking. The pre-heating 
has little or no effect on the mass of gas that cools until the pre-heat temperature reaches 10 7 ' 9 K, 
about 25% of the perturbation's virial temperature. At higher temperatures, pressure support 
prevents gas shells from turning around and collapsing. 

Allowing gas a Hubble time to cool may be unfairly generous for a hierarchical scenario. A 
perturbation will often merge with another of comparable size only a few dynamical times after it 
collapses, and gas that has not cooled by then can be shock-heated to a higher temperature. The 
appropriate time to allow for cooling may therefore be a few dynamical times rather than a Hubble 
time. The triangles in Figure 12 show, for several pure fluid calculations, the gas that cools within 
a single dynamical time, i.e. a gas shell that shocks at t s must cool by time t s + (37T/16/0) 1 / 2 , 
where p is the average interior density immediately after the shock. Since there is less time for 
cooling (and the ordinate M s is still the mass that shocks by z = 0), the cooled masses are lower. 
However, while the trend of M c vs. M s is somewhat shallower, it is certainly not flat. Even the 
stringent requirement of cooling within a single dynamical time cannot by itself produce a sharp 
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cutoff in the galaxy luminosity function. 

Radiative cooling is bound to be an important ingredient in the process of galaxy formation. 
Gas must cool before it can form stars, and dissipation is needed to explain the prevalence 
of galaxy disks, to explain the difference between the characteristic mass of galaxies and the 
characteristic mass of rich clusters, and, in most scenarios, to explain the high internal densities of 
ellipticals and bulges. However, our results indicate that the characteristic mass of galaxies cannot 
simply be read out of the physical constants that describe gravity and atomic physics. Instead, the 
form of the galaxy luminosity function must reflect a more subtle interplay between cooling and 
cosmology, particularly the rates at which perturbations collapse and merge. While this interplay 
complicates our effort to understand the physics of galaxy formation, it raises the hope that the 
properties of galaxies may provide important constraints on cosmological parameters, the nature 
of dark matter, and the spectrum of fluctuations in the early universe. 
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